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We discuss the high density behavior of a system of hard spheres of diameter d on the hypercubic 
lattice of dimension n, in the hmit n 00, d — > 00, d/n = 5. The problem is relevant for coding 
theory, and the best available bounds state that the maximum density of the system falls in the 
interval 1 < pVd < exp(nK(5)), being k{S) > and Vd the volume of a sphere of radius d. We find a 
solution of the equations describing the liquid up to an exponentially large value of p = pVd , but we 
show that this solution gives a negative entropy for the liquid phase for n. We then conjecture 
that a phase transition towards a different phase might take place, and we discuss possible scenarios 
for this transition. 

A code C is a subset of the binary Hamming space H2 = {0, l}". The distance between two points s,t G i?J is 
the Hamming distance d{s,t) — + ia)mod2, i-e. it is given by the number of different bits. We consider the 

problem of finding the maximal size A(n, d) of a code C such that the minimum distance between two points in C is 
d, that means, denoting by \C\ the number of sequences in C, 

= max [|C| Ws,t & C, d{s,t) > d] . (1) 

In particular we are interested in the quantity 

R{6) = limsup - \og2 A{n, d) , (2) 

n — >oc, d/n — y6 ^ 

where the supremum is taken on all possible sequences of codes such that d/n — s- 5. The problem trivializes for 
5 > 1/2 as the total number of sequences is finite and R{5) = 0. An interesting scaling is d/n = 1/2 — en~", as for an 
appropriate choice of a the number A{n, d) might increase polynomially in n, but this will not be investigated here. 
Thus we will restrict to J < 1/2 in the following. 

This problem is relevant for the theory of error correcting codes [1-6] . In "physics language" , it is the problem of 
finding the maximum possible density of a system of hard spheres on the hypercubic lattice. This rephrasing of the 
problem has been shown to be useful as it allows to use well known methods borrowed from the theory of liquids, like 
the virial expansion [4]. 

In this paper we will discuss the behavior of the system at high density, in order to understand how one can try 
to compute the maximum density. We will show that, for large n, the problem closely resembles the problem of hard 
spheres in M" in the limit of large space dimension n. The basic idea is that in this limit the number of neighbors of a 
sphere is large, much as it happens in the continuum for large space dimension. We will then discuss some recent ideas 
that have been used in the continuum [7-10] to make some progress in the direction of deriving bounds on A{n, d). 

The best known lower bound on A{n,d) (Varshamov- Gilbert bound) states that the density p = |C|/2" > 1/Vd-i, 
Vd being the volume of a sphere of radius d in H2 [1] . This bound can be proven from the convergence of the virial 
series [4] and gives R{S) > Rvg{S) = 1 — H{5), H{5) being the binary entropy function (see below). This means that 
a "liquid phase", defined by the virial equation of state, exists at least up to pVd-i ^ 1. We will show that the hquid 
phase can be formally continued up to a density pVd~i ^ exp{nK{S)) with k{S) > 0: this result correspond to the 
Frisch-Percus result in the continuum [7] . 

However, we find that the entropy of the liquid becomes negative at pVd-i ^ n. This suggests the possible instability 
of the liquid towards a different phase, i.e. the existence of a phase transition. We will discuss two different possibilites. 
By analogy with the problem in R for large n, we conjecture that a glass transition might be present also in this 
system. In absence of other phases (such as "crystalline" phases), this analogy suggests that R{S) is given by the 
Varshamov-Gilbert result, R{6) = Rvg{S)- We will also discuss a different instability that happens for even d, leading 
to a first order transition to a phase where the particles move on a sublattice of H2. 

The calculation of the properties of the (eventual) glassy phase requires the use of the replica method, but this 
turns out to be more difficult than in M so we leave it for future work. 

The paper is organized as follows: in section I we set up the basic notations and definitions; in section II we review 
the known bounds on R{S)] in section III we present our main results and conjectures on the high density behavior 
of the system; finally in section IV we draw the conclusions and present a summary of our ideas. 
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I. HARD SPHERES ON THE HYPERCUBE 
A. Definitions 

Let s = (si, • • • , Sn), Sa & {0, 1}, be a point in H2. We define the Hamming distance between two points 

n 

d{s, t) = '^{Sa+ ia)mod2 • (3) 
a=l 

The number of points at distance i from the origin is Si = (") and the volume of a sphere of radius d centered in the 
origin is 

d 

' n 
i 



t:d(0,t)<d i=0 



The interaction potential between two particles in s and t is the hard core potential: 



U{s,t) = i (5) 
I oo a(s, t) < d . 



The grancanonical partition function [4, 11] is 



2" jsi 
Z 



^"W=Em e-^.<.-(-^), (6) 

the average number of spheres (to be identified with the average size of the code, |C|) is (N) = and A{n, d) = 

lim2_>oo {N). The density is p = (A'') /2" and it is convenient to define the reduced density 



{N} Vd. 

and 



p = pVd-1 = — , (7) 



¥'=^lnp. (8) 



For n ^ oo and d/n = 6 < 1/2 one has 



v{6) = lim - InVns = -6ln6-{l-5) ln(l - 5) , (9) 



i.e. yd~.5d = Q ~e™W, so that 

y,= llnp~ ^lnp + t;((5) , p ~ e"!*'-''^'')! . (10) 

The pair distribution function is given by the average number of particle pairs such that one particle is in s and the 
other in t [11]: 

P2{s,t) = /^Ss,,,St,s,) = p^9{d{s,t)) (11) 

where the Kronecker 5s,t is equal to one iis = t and to zero otherwise. The pair distribution function g{i) is normalized 
to have value 1 at large distance (for n — > oo). We also define h{i) = g{i) — 1 which vanishes at large distance. Note 
that the function pg{d{0, s)) is the probability of finding a particle in the point s, given that there is a particle in the 
origin. 
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B. Liquid phase 



The liquid phase is defined by the virial expansion of the partition function at low densities which converges for 

pVd-i < l/6e [4]. However, one can try to perform a contimiation of the liquid equation of state below the radius 
of convergence of the virial series. For hard spheres in , such a continuation is possible: approximated expressions 
are usually obtained by resumming some class of diagrams of the virial expansion to get a closed integral equation 
for the pair correlation function g{i). Well known resummations are the Percus-Yevick (PY) and the HyperNetted 
Chain (HNC) ones, and agree well with numerical results [11]. 

For hard spheres in R" it has been shown by Frisch and Percus [7] that in the limit of large n the virial series is 
dominated order by order by the so-called ring diagrams. The resiimmation of the ring diagrams has been shown to 
provide a reasonable analytic continuation of the liquid equation of state up to very large values of the density [7] . As 
the latter diagrams are included in the HNC resummation, the two resummations should be equivalent in the large 
dimension limit. The advantage of the HNC resummation is that it leads to a closed expression for the free energy 
corresponding to the partition function (6) for any finite n. The HNC free energy (per particle) is: 



m = f E Ui) ln5(i) - 9{i) + 1] + - 1 + ^ E ' (12) 

i=0 ^ p>3 ^ 

where 

TihF = E Ksi, S2)h{s2, si) ■ ■ ■ h{sp, si) , (13) 

Sl •••Sp 

and the function g{i) is determined by the stationarity condition ^ = {HNC equation). We will argue that a 
solution to the HNC equations (or to similar approximations to the equation of state of the liquid, that should be 
equivalent for n — > oo) exists up to p ~ exp(nK) with k > 0, thus defining the continuation of the liquid equation of 
state in this region of densities. 

C. Fourier transform 

To write the last term of the HNC free energy in a more convenient way it is useful to define the Fourier transform 

on the hypercube. Define the scalar product of two sequences as s ■ t = X]q=i ^ata (such that s"^ = s ■ s = d{0, s)) 
and let us indicate by s + t the sum (modulo 2) of the two sequences. Then if q G {0, 1}" the Fourier transform of a 
function h{s) is given by 

%) = E(-i)'"'^(^)' 

(14) 

2^' 



q 



Using the property X)j,(— l)*-^^^''^-''* = 2'^5q^_q^ the Fourier transform of a function /i(s, t) = h{s + t) is 

^(_l)«i-(_l).-t/,(5 + t) =2%,,,,%i) , 

s.t 

(15) 

2" 



h{s,t) = ^Y.^-ir^^+'^'^^q). 



Moreover if h{s,t) = h{d{s,t)) (a rotationally invariant function), its Fourier transform depends only on a = q^ 
d{0,q), see Eq. (3), and one has 



h{a) = '^J^n{a,i)h{i) , 

1 " 



2^^ „ 

a=0 

min(a,2) 

Ma,i)= E (-1)^--= E i-^r (:)(:::)= 

s:d{0,s)=i m— max(0,i+a— n) 



(16) 
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where Ki{n,a) is a Krawtchouk polynomial [6, 12]. The matrix ^„(a, i) has the following symmetries that will be 
useful in the following: 

□ .F„(a,*) = ('^).F„(*,a) , (17) 
J^n{n- a,i) ^ {-'iyj^n{a,i) , (18) 
J^n{a,n-i) = {-irj^nia,i) , (19) 

-a,n-i) = (-l)"-«-^jr„(a, i) . (20) 

It satisfies ^„(0, i) = (") and JF„(a, 0) = 1, and can be easily constructed using the recursion relations 

J'„+i(a,0) = 1 , 

J^n+ii.a,i + i) = ^n{a,i) +^n{a,i+^) , a<n, (21) 
:F„+i(n+ l,i + 1) = -^„(n,i) +:F„(n,i + 1) . 

Using Eq. (15) it is easy to show that 

Trh^ = Y^[hiq)r , (22) 



and the HNC free energy (12) becomes, defining L3{x) = ln(l + x) — x ■ 



2 ' 



n _j n 

= f E ^^9(i) - 9{i) + 1] + Inp - 1 - j— ^ S,Ls[ph{a)] , (23) 

i=0 ^ a=0 

and the HNC equation can be written as [11] 



ln5(s) = h{s) - c{s) , 
1 + ph{q) 



II. KNOWN BOUNDS ON R{S) 

In this section we will review some known bounds on R{6). Using Eq.s (2), (10), and p = (N) /2", R{5) is related 
to the maximum density <fc{S) of the spheres by 

R(5) = l-H{5) + ^ , (25) 

where H{6) = = — Jlog2 6 — {I — S) log2(l — 6) is the binary entropy function. We will use units of ip because it 
will lighten the notation in section III. 

As discussed in the introduction, the best lower bound for R{6) (Varshamov-Gilbert bound [1]) can be obtained by 
proving the convergence of the virial series [4]. It turns out that the virial series converges for pVd-i = p < l/6e, that 
means </? < for n — > oo. Thus (fc>0 and a lower bound for R is R{6) > 1 — H{6). 

A trivial upper bound follows from the fact that the total volume occupied by the spheres is ^ NV^i2 and should 
obviously be smaller than the total volume 2". Thus 

= S < t7- ~ ^""'^'^'^ ^ 'fi<Vma. = v{5)-v{5/2) , (26) 
2" Vd/2 

so 1 — H{5) + ^max/ ln2 = 1 — H{6/2) is an upper bound for R{5). 

Better upper bounds for R{8) can be derived by Delsarte's linear programming method [2]. In physics language, 
they follow from the observation that the minimal requirements for the correlation function g{i) are the following: 

i) g{i) = for < i < d, as no pair of particles can be at a distance smaller than d due to the hard core interaction. 

ii) 9{i) > 0, Vi ; this follows from the definition of g{i), see Eq. (11). 

Hi) ph{a) > —1, Va; this is because the structure factor S{a) = 1 + ph{a) is a positive quantity, equal to the average 
of the square modulus of the q-component of density fluctuations [11]. 
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FIG. 1: Bounds for tpc- The VG lower bound is > 0. The upper bounds v'™ax(5) (26), dashed hue, and ifMRRwi^) 
(32), full line, are reported. The lower bound to D{5) obtained by Samorodnitsky, vps(5) = 'Pmrrw{S)/2, is reported as a 
dot-dashed line. Finally, the value of ipHNc{S), below which our solution for the liquid state is well defined, is reported as a 
dot-dashed-dashed line. 



If one is able to find a function 17(1) verifying conditions above (that have sometimes been called positivity 

conditions, see e.g. [10]) for a given value of p, one can compute the corresponding value of (N) either as {N) = p2" 
or by working in the canonical ensemble (i.e. at fixed N [11]) and recalling that from Eq. (11) it follows: 



N{N -1) 



s.t 



1=0 



P9{i) 



N = 1 + 



E 

1=0 



(27) 



It is not obvious that one can find configurations of the system with density p that actually produce the function 
g(i), see e.g. [10, 13, 14] for a discussion of this issue in the case of spheres in the continuum. However, as the p, g(i) 
obtained from the partition function (6) (or from the canonical partition function) must satisfy conditions i)-iii), it 
is clear that the value of A{n, d) is smaller than the maximum of the right hand side of Eq. (27) over all the possible 
choices of p, g{i) satisfying the positivity conditions, i.e. 



A{n, d) < max 



1 



71 

E 

1=0 



(28) 



The problem with Delsarte's method is that it can be used to derive arithmetically an upper bound for finite and not 
too large n (e.g. n = 1000, see [5]), but it is not easy to obtain analytical results for n 00. 

In the literature the Delsarte method has been often formulated in terms of the function Ai = Sio + p(J^^g{i), see 
e.g. [5, 6]. Conditions i)-ii) are equivalent to ^0 = 1, = for 1 < i < d and Ai > 0. Condition Hi) gives, using 
Eq.s (16) and (17), and recalling that J>i(a, 0) = 1, ^„(0,i) = (") and X]"=o •^"('^' *) ~ 2"(5ao (being the Fourier 
transform of the function 1): 



5(a) = l + ph{a) ^1 + J2^nia,i)[pg{i)- p] = 1 + E 



^Tn{i,a) 



i=0 



p\ ■ ]gii) + So'i. 



prSao 



J^nji, a ) ( n 

n 
a 



pg{i) - pT'Sao 
Y,^n{i,a)A,- p2'^5ao ■ 



(29) 



Thus the condition S{a) > is equivalent to X]"=o -^"(^j '^)^« > for a ^ 0, while for a = we simply have 
Sr=o = from Eq. (27). Then the bound (28) can be reformulated as 



A{n, d) < max 



i=0 



'A, 



A, > Vi ; Ao = 1 ; Ai 



i=0 



D{n,d) 



(30) 
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Upper bounds for A(n.d) can be derived by studying the dual to the linear problem (30), see e.g. [5]. In this way 
one can prove that (MRRW bound [3]), for n ^ oo: 

D{S) = hmsup - log2 D{n, d) < H - ^/6{l - 6)] = Rmrrw{S) , (31) 

that means that R{S) verifies the same bound and 

^ <y (^l - ,/S(]rrs)^ -ln2 + v{S)=ipMRRw{S) . (32) 

This bound has been (little) improved only for S < 0.273 [3, 5]. 

By constructing an explicit solution to the positivity conditions, one can obtain lower bounds on D{S) (but this is 
not sufficient to have the same bound for R{d), [10, 14]). For this problem, this has been done in [6], where it was 
shown that 

g{i) = {d+l)5id + 0ii-d) , (33) 

where 9{i — d) = 1 for i > d and otherwise, verifies conditions up to a density ifis{S) = <Pmrrw{S)/2. 

To resume, the situation is the following: 

• the VG bound states that the maximum density (^c > 0; 

• the MRRW bounds coming from the dual to (30) prove that (pc < fMRRw{S), Eq. (32) (with a little improvement 
for S < 0.273, see [5]); 

• the Delsarte's method cannot be used to prove that the VG bound is tight, i.e. that ipc < 0, as the lower 
bound on D{S) obtained in [6] implies that the best upper bound from Delsarte's method cannot be smaller 

than (pMRRw{S)/2; 

• moreover, recent arithmetical results for n = 1000 [5] seems to indicate that the actual value of D(5) might be 
close to the MRRW ones, i.e. much larger then the lower bound of [6]. This means that the MRRW bounds are 
probably very close to the best one can obtain from the Delsarte's method (30). 

The bounds above for ipc (summarized in Fig. 1) leave a large gap - at least of the order of ^pmrrw{5) /2, but probably 
of the order of <Pmrrw{5) - and there are not so many ideas on how to improve them [5]. Note that, as discussed in 
the introduction, for 5 > 1/2 one can prove that R{5) = 0, so the size of the code is not exponential. 



III. THE LIQUID PHASE AT HIGH DENSITY 

In this section we will discuss some insight on the problem that comes from the physical intuition on the possible 
behavior of the system (6). We are not able to present rigorous results but we hope that the discussion below will 
lead to new ideas on how to rigorously improve the bounds on R{S). 

It is convenient to outline our basic ideas before going into the details of the calculations. We try to find a solution 
to the HNC equations (12), (24) (or to other approximate equations for the liquid) for > {(p is defined in Eq. (8)). 
We assume here that any resummation will be equivalent for n — » cxd as long as it includes the ring diagrams [7]. 
Such a solution should clearly verify at least the positivity conditions However there can be many different 

solutions to these conditions that may not correspond to the high-density liquid. In particular the solution proposed 
by Samorodnitsky, Eq. (33), is not suitable to describe a liquid state, as we do not expect to observe a large number 
of particles in contact (represented by a peak at i = d) in the liquid phase. Thus we will first look for a function g{i) 
verifying i)-iii), not showing large peaks and departing continuously from the step function for > 0. We will show 
in the following that such a function exists up to = ifiHNciS) < fsi^), see Fig. 1, and is indeed given by the step 
function plus an exponentially small correction in n. We interpret this solution as describing the liquid phase and 
show numerically that the solution of the HNG equations converges to this solution for n — > oo. 

Using the solution above we can compute the entropy of the liquid. A crucial observation is that this entropy 
becomes negative for ip ^ \og{n)/n, i.e. very close to the VG bound. This means that the liquid phase must become 
unstable below this value of density, as the entropy of a discrete system must be positive. We then expect that the 
system (6) will undergo a phase transition at a density ip < log(n)/n. 

This behavior closely resembles the behavior of hard spheres in the continuum in the limit of large space dimension. 
For this problem, we recently showed [9] that at a value of density close to the radius of convergence of the virial 
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expansion (i.e. to the VG value) the liquid phase becomes unstable towards a glass phase where replica symmetry is 
broken. In the glass phase the pressure rapidly increases and diverges at a maximum density (for the glass) which is 
found to be of the same order of the glass transition density. 

By analogy with the problem in the continuum, we argue that also in this problem a glass transition exists at a 
density (pK < log(n)/n and that the glass phase should exist up to a maximum density with the same scaling in 
n. This means that, if no other phases exist, one should have ifc ~ log(n)/n, i.e. the VG boimd should be tight. 
Unfortunately we are still not able to repeat the calculation of the equation of state of the glass - that was done in 
[9] for the problem in the continuum - for this problem; but we believe that the computation is feasible and we leave 
it for future work. 

Clearly other phases may exist at least for some special values of n, d. For instance, we have some numerical 
evidence, for even d, of a first-order phase transition toward a phase in which particles only occupy a subspace of the 
Hamming space. We will discuss this issue below. 

In the following we will try to make these arguments more precise. 



A. The Fourier transform for n — » oo 



We begin by studying the properties of the Fourier transform of the delta function g{j) = Sij, that is simply ^„ (a, i), 
for n ^ 00. We define k = a/n and x = i/n. We want to compute the function f{k, x) defined by 



<^{k,x) = J^n{nk,nx) ~ Ree"^^^'^) , n ^ oo . 



(34) 



Using the fact that the poles of the gamma function T{z) are in ^ = —k with residual (— l)'^/^!, we can rewrite 
Eq. (16) as 



J^n{a,i) = / dzT{z) — 
Jc ^ \C' 



r(a + l)r(n - a + 1) 



\a + z + l)T{i + Z + l)r(n -a-i-z + 1) ' 



(35) 



where the contour C embraces the negative part of the real z axis, see Fig. 2. Using the Stirling formula r(n) 
n"e~"\/27rn, we get, changing the integration variable to C = z/n, and neglecting power-law prefactors. 



Tn{nk,nx) = / d^e 
Jc 



,n[C logC-tj(/c)-(/c+C) log(fe+C)-(a;+C) log(x+0-(l-k-x-C) log(l-fe-a:-C)] 



/ dCe""^'^'^ 
Jc 



and we can evaluate the integral using the saddle point method. The saddle point equation is 

a^-k-x-0 



{k + c){x + o 



1 



with solutions 



--k-x 



±-\ k-- 



X- - - - 



(36) 



(37) 



(38) 



From the analysis of the position of the solutions in the complex plane one can deduce the following: 



1. In the region {k — + {x — — j < (region A in Fig. 2) the solutions are complex, and C+ = C-> 

■ /(fc,x) = a(C-) . (39) 



cr(C+) = [cr(C-)]*. Thus One has 

J^n{nk,nx) ~ e"'^^^-) + e""(«+) = Ree"^^^-) ~ enRea(C-) cos[nIma(C-)] 



2. In region B of Fig. 2 the saddle points are real and positive with < C- < C+i and cr(C±) are also real with 
cr(C+) > fT(C-) > 0. The point (- is then the closest saddle point to the original integration contour; moreover, 
it is a local mimimum of a{() along the real axis, so it is a maximum of a{() along the imaginary direction and 
the integration path can be deformed to include it without crossing regions of C where Recr(C) > o'(C-), see 
Fig. 2. Then 



J^„{nk, nx) ~ e"''^^-) ^ f{k, x) = a(C-) 



(40) 
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FIG. 2: Left: the domain in tlie plane (fc, x) where the saddle points (38) are complex is represented by region A. Right: 
deformation of the integration contour to include the saddle point C- (this plot refers to region B, see text). 



3. The behavior in the regions C,D,E can be obtained using the symmetries (18), (19), (20). Alternatively one can 
always choose the closest saddle point to the integration path that is a local minimum on the real axis: it turns 
out that one has to choose C- in the region E and C+ in the regions C and D. With this choice the symmetries 
(18), (19), (20) are respected. 

Finally one obtains 

f{k,x) = mm[a{C+),a{C-)] , (41) 

Re 

where min^e means that one has to take the solution with the smallest real part. The real part of the resulting 
function f{k,x) is an increasing function of x for all k and x < 1/2, see Fig. 3. This allows to compute the Fourier 
transform of the theta function 9{d — i) for d < n/2. Defining /a;(fc, x) = dxf{k, x), we have 



es{nk) = Y^Tn{nk,i)e{d - i) nRe / rfa;e"^('='^) ~ nRe / rfa;e" 

,_n ^0 Jo 



[f{k,6)+f^{k,S){x-&)+-] 



i=o - - (42) 

e"/(fe.5) 

i.e. the Fourier transform of the theta coincides with the one of the delta to leading order in n as long as d < n/2. 
Similarly we can compute the Fourier transform of a function that vanishes outside a finite interval and approaches 
zero linearly at the edge of the interval, i.e. h{i) = 0{d — i){d — i): 

h{nk) r-. n'^Re f dxe"f^'''''\S - x) n^Re [\xe"^f^'''^'^+f''^'''^^^^-^'>+-'^{5 - x) 

Jo Jo (43) 

i.e. to leading order in n also this function is equal to $(&, S). 



B. HNC for n -» oo 



We argue that for n ^ oo the solution to the HNC equation approaches a solution of conditions i)-iii) which shows 
no large peaks. In particular we will look for a solution of the form g{i) = 6{i — d + + exp{nh{i))], h{i) < 0, i.e. 
a solution differing from the step function (that describes the liquid for < 0) by an exponentially small quantity. 
Conditions are, for x = i/n and k = a/n, and recalling that h{i) = g(i) — 1, 

' ph{k) > -1 , 

< h{x) > -1 , (44) 
h{x) = — 1 , X < S . 
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FIG. 3: The function Ref{k, x) for different values of fc as a function of x; from top to bottom k ■ 
that from Eq. (18) Re/(A;,a;) = Re/(1 - A;, a;). 



0.1,0.2,0.3,0.4,0.5. Recall 



First we will check that for < the step function g{x) = 6{x — 5), i.e. h{x) = —9{S — x), satisfies the conditions 
above. Using Eq.s (10) and (42), 

ph{k) = -p$(fc,(5) > -1 ^ ^ee''^'P-viS)+f{k,S)] ^ ^n[v-vi5)+Ref(k,S)] cos[nlm f {k , 6)] < 1 , (45) 
and the latter relation is equivalent to 



Rc fik,d) < v{S) - ip 



(46) 



In Fig. 4 the function Re f{k, 6) is reported as a function of k for a representative value of S. It assumes its maximum 
in fc = and k = 1 and f{0,6) = v{6). Thus the inequality (46) is always satisfied if (p < 0, so that in this region 
(which is also the region where the virial scries converges) wc argue that h(x) = —9(5 — x) describes the liquid phase 
for n oo. This can also be checked by a direct evaluation of the leading terms in the (convergent) virial series, 
see [7]. Following [7] we also argue that the HNC resummation contains all the relevant diagrams for ip < 0, so we can 
use it to obtain the free energy of the liquid. Substituting the result for h{x) in Eq. (23) the last term is exponentially 
small in n and one obtains, up to exponentially small corrections. 



p dp 2 



(47) 



where P is the reduced pressure. As found in [7] wc find that the entropy is given by the ideal gas term plus the first 
virial correction. Note that p is exponentially small for p < so the system behaves essentially as an ideal gas. 

For ip > the function h{x) cannot be given by —9{d — x) as this function does not respect the positivity conditions 
(44). We follow the strategy of [8] and decompose h{x) = ho{x) + hi{x) assuming that ho{x) vanishes for a; > i5 and 
hi{x) is a continuous function of x. We call Y = hi{S) = — 1 — ho{S~), i.e. ho{S~) = — (1 + Y). As ho{x) vanishes for 
x > 5, its Fourier transform is given by (42): 



ho{k) = -{l + Y)^k,6) , 
and the condition ph{k) = pho{k) + phi{k) > — 1 becomes 

hi{k) > -- + {1 + Y)<b{k,5) 
P 

We choose the simplest solution to the previous equation, 



for - i + (1 + r)$(/c, (5) > 
for - i + (l + r)$(A;,(5) < 



(48) 



(49) 



(50) 
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FIG. 4: The function Re/(/c, S) as a function of for J = 0.1. The value of v{S) and of v{6) — 0.05 are also reported. The value 
of k such that Re f{k, 6) intersects v{5) — (p is the value of A;c(i5, ip), see Eq. (53). 

and wc will show that it gives, in real space, an exponentially small correction with respect to the step function. The 
function hi{k) has a two singularities when 

(1 + F)p$(fc,(5) = 1 ^ ip-vi6) + fik,S) + -ln{l + Y)=0, (51) 

n 

but the above equation is well defined only if its solutions lie in the region where /(fc, 6) is real, otherwise the function 
$(fc, 6) will oscillate very fast. This will impose some restrictions to the values of ip and d. 

First we will restrict to odd d in order to have the symmetry ^{1 — k,S) = (f>(fc, (5): for even d we have the opposite 
symmetry $(1 — k,5) — — $(fc, S) and f{k, 5) has an imaginary part for k > 1/2. This follows from <I>(fc, 6) — !Fn{a, d) 
(with a = nk and d = nS) and from Eq. (18). For odd d Eq. (51) will have two solutions kc, 1 — fcc due to the 
symmetry, see Fig. 4. Note that the opposite restriction was applied in the numerical computation of [5] where only 
the case of even d has been considered. 

Next, we look for a solution k^. outside the region A of Fig. 2, as in region A we already know that f{k, S) is not 
real. This means that the maximum possible value for kc is 

fcma.(<5) = ^-y'i- (.5-0 =1-7^(13^, (52) 

that is the boundary of the region where /(fc, S) is real, see Fig. 2. 

As we will self-consistently verify at the end, under the restrictions above one has Y ~ e""", then Eq. (51) becomes 

f{k,6) = v{6)-ip . (53) 

For < it has no solutions as discussed above, so hi{k) = 0, then hi{x) = 0, F = and we recover the step function 
solution. For ip > 0, kc increases from to kmax{S) and reaches the boundary of region A, using Eq. (53), exactly at 

ip = v{6) - f{kmax{S), 5) = ipMRRw{5)/2 = ips{5) , (54) 

where ^Pmrrw{^) has been defined in Eq. (32). 

We have now to compute hi{x) from hi{k). The function hi{k) verifies hi{k) = hi{l — k) and is nonzero only for 

k < kc and k > 1 — kc, it vanishes linearly close to kc, hi{k) ~ (1 + Y)^k{kc, S){k — kc), and ^k{k, S) = dk^{k, 5) is 
real in fc = kc- If we restrict to even i, we have also ^{x, k) = $(x, 1 — A;) by symmetry (19). We can do that because 
it is possible to show that, in the case of odd d that we are considering, one has h(i — 1) = h(i) for even i; thus it is 
enough to compute h{i) for even i (see the next section, Appendix A and Fig. 6 for a detailed discussion of this tricky 
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point). Using Eq. (43) we have then 



2r? /"^/^ 2n P^" 

^i(^) = ^y^ dA:$(a;,fc)/ii(A:)~ — (l + y)e"^('=-^)n/fc(/ec,5)Re d/ce"/(^>'=)(A: - fee) 



(55) 



2"-i [/.(cc,fcc)]2 ■ 

Keeping only the leading terms (exponentials of n) we obtain the self consistency equation for Y: 

Y = hi{S) = (1 + y)e"[/(fcc,5)+/(5,fc.)-in2] ^5g^ 

If fee < fcmaa; wc have /(/cc, (5) + /((5, fc^) — In 2 < 0, SO that y is exponentially small and is given by 

Y = gn[/(fec,<S) + /(<5,fce)-ln2] ^ ^n[v{S)+vikc)-2^-in2] _ ^gr,^ 

where we used the relation f{x,k) = ,f(k,x) + v{k) — v{x) that follows from Eq. (17) and Eq. (53). Finally, the 
function h{k) = ho{k) + hi{k) is given, using Eq. (50) and F ^ 1, by 

- J -1 for fc < kc and k > 1 - kc , 

We can rewrite Eq. (55) using F <C 1, neglecting non-exponential prefactors, and recalling that f{kc,5) = v{5) — ip, 
as 

^n[v{S)-ip] I I 

hi{x) = — ^{x, kc) = -^^(a;, kc) ; (59) 

note that from Eq. (16) it follows that 2~"$(a;, kc) is the Fourier transform of 5{k — kc). Finally, h{x) is given, from 
Eq. (55), by 



h{x) 



-1 X <S , 

^^Q(.nHx) x>5, (60) 

h{x) = v{6) -ip + f{x, kc) - In 2 . 



The solution above is defined up to the value ipg of the density given by Eq. (54). Indeed, the solution fc^ of Eq. (51) 
is given by Eq. (53) only if it is in the region where f{k,6) is real. Otherwise, oscillations are present in $(A:,x) and 
the solution is not well defined. Moreover, if Y is not exponentially small again the solution above fails. Both these 
conditions seem to be violated for 1,5 > 1^2.5 ((5). 

There is however another condition to be imposed, namely that h{x) < 0; otherwise the solution will be exponentially 
large and again we do not expect that for a liquid phase. The maximum of h{x) is attained in a; = 0, x = 1, and is 
given by v{6) + v{kc) — — In 2. The condition h{x) < then requires 

v{6) + v{kc) - V - In 2 < . (61) 

A numerical solution of the previous equation (recall that kc depends implicitly on <f) gives the stability threshold 
(pHNciS)^ which is reported in Fig. 1. For (p > (fHNC our solution starts to exhibit diverging oscillations at large x. 
This is observed also in the numerical solution of the HNC equations, see below. Above <^s((5) either the solution does 
not exist anymore or it yields a value of Y that is exponentially diverging with n. In both cases the solution does not 
describe a liquid phase. We will see however that we are not really interested in so high values of the density as the 
liquid phase becomes unstable at a much lower density. 



C. HNC entropy 

The HNC free energy is the canonical free energy, that for hard spheres is simply —S{p). For < it is given by 
Eq. (47), and we showed that for > only exponentially small corrections appear. Neglecting these corrections we 
have from Eq. (47) 

S{p) = 1 - logp - ^pVd ~ 1 - + nv{S) - ^ log[2^(5(l -S)]-^ logn - ^e"'^ , (62) 
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FIG. 5: Comparison of the numerical solution of the HNC equations with the asymptotic solution of the positivity conditions. 
Left: The function \n[\ph(k)\]/n obtained numerically is reported as a function of k for n = 36 and d = 11 (i.e. 5 — 0.305) and 
if = 0,0.05,0.1. The lines are obtained from Eq. (65). Right: The function h{i) for n = 36, d = 11 and ip = 0,0.05,0.1,0.12. 
At large (p we start to observe diverging oscillations at high values of i, see text. 



in the full range if < fHNci^)- An interesting observation is that S{p) becomes negative at a value of density given 
by, keeping only the leading terms, 

S{ip) - -nip + nv{S) - i logn - ^e"-^ = , (63) 
and the solution is, at first order in n oo, 

Lpo = - log[2nv{S) - 3 log n] . (64) 
n 

As for a discrete system S{ip) > 0, this means that the liquid phase must become unstable at a density ip < ipo ^ 
log(n)/n. 



D. Numerical solution of the HNC equations 

We will now compare the asymptotic solution with a numerical solution of the HNC equations (24). The latter are 
solved using a standard iterative algorithm. For graphical convenience we report the function ^ log \ ph{k)\, that from 
Eq. (58) is given by 

r li \h(k]\-}'^ ioT k < kc and k > 1 ~ kc , , . 

n'^n ^ \Ref{k,S)+^-v{S) + ^\og\cosnlmf{k,S)\ ioi k, < k < 1 - k, . ^ ^ 

The last term is nonzero only for k S [kmaxi^), 1 — kmax{^)] and gives rise to oscillations whose frequency increases on 
increasing n. Moreover at the values of k where cos nlm f (k , S) = the last term diverges. These values accumulate 
in the interval k g [kmax{S), 1 ~ k„iax{^)] for large n. However, for the relativaly small values of n considered here, 
the are no integer values oi a = nk where cosnIm/(/c, (5) = and we can compare the function ^ log \ph{k)\ with its 
asymptotic limit neglecting the last term. Moreover we can compare h{i) with the analytical expression (60). The 
comparisons are encouraging as shown in Fig. 5. We find that h{k) agrees well with Eq. (58) and, for small ip, h{i) is 
essentially the step function plus a small corrrection. On increasing the density large oscillations appear at large i, as 
predicted by Eq. (60). Unfortunately a quantitative comparison of h(i) with the asymptotic expression requires either 
the evaluation of finite n corrections, due to the small values of n we can investigate numerically, or the (difficult) 
investigation of much larger values of n, see e.g. [5]. 
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FIG. 6: Left: correlation function h{i) for n — 36, d = 10 and ip — 0,0.016. Inset; structure factor S{a) = 1 + ph(a) for the 
same values of the parameteres. On increasing above = strong oscillations appear in h{i), related to a strong peak in 
S{a) for a — 36. Right: same plot for n = 36, d = 11 and <p = 0,0.05. In this case the oscillations are absent, S{a) has a 
different behavior for large a, and h{i) decays in steps of two. The different (approximate) symmetry of S{a) is related to the 
symmetry (18). 

E. Is there a glass transition? 

The fact that the entropy becomes negative seems to indicate the existence of a phase transition. By analogy with 
the continuum problem, where we showed that a glass transition happens at a similar value of the density [9], we can 
conjecture that a glass transition will happen also in this problem. 

If this is the case, one can show by general arguments and by analogy with the continuum problem that S glass (v) ^ 
Siiquidif) (this is because a downward jump of the compressibility is expected at the glass transition [15]). This 
means that the entropy of the glass will vanish at a density (p < ipo \og{n)/n. Then both the liquid and the glass 
phase will disappear for (p ^ log(n)/n. If there are no other phases (such as a "crystal"), this scenario indicates that 
the Varshamov- Gilbert bound should give asymptotically the exact result. 

Note also that the correlation function used by Samorodnitsky, Eq. (33) is similar to the one we expect for the 
glass phase (see the discussion in [9, 15]). This means that, if the picture above is correct, the g{i) in Eq. (33) should 
correspond to realizable packings only up to a density (f ~ log(n)/n. We hope that further work will clarify this issue, 
see also the discussion in [10, 14]. 

F. Other instabiUties 

We would like now to consider a different instability of the liquid that is observed for even d. Indeed, d is the smallest 
possible distance between any pair of spheres. Let us split the Hamming space i/^ into two spaces H2 = -ff" U i/", 
where 7?^" = {s € H^\d{s,0) is even} and = {s € H^\d{s,0) is odd}. Recah that pg{i) = p[l + h{i)] is the 
probability of finding a particle in a given point at distance i from the origin, given that there is a particle in the 
origin. 

At high density, given that there is a sphere sq in the origin, it will be convenient, say, to place the second sphere 
si at the minimum distance d, thus si € H^. Then we can put another sphere S2 at minimum distance d from si, 
so that S2 G H^, and so on. It is clear that this picture is oversimplified but still we can expect that at large enough 
density particles will concentrate on one of the subsets H^^, depending whether there is or not a particle in the origin. 

This is indeed what we found in the solution of the HNC equations. The partial localization of the particles on a 
subset -ff"o is revealed by oscillations in h(i), which is much bigger on the even values of i, see Fig. 6. The oscillations 

in h(i) are related to a strong peak developping in S{a) = 1 + ph{a) for a = rt on increasing tp above 0. Due to 
this growing peak, the HNC equations rapidly become unstable for small values of p (much smaller than for odd 
d). Preliminary Montecarlo simulations seems to indicate the existence of a first order phase transition to a phase in 
which particles are completely localized on a sublattice at a critical value of (p. 

For odd d a related phenomenon occurs. Indeed, by looking at the correlation function h(i) for odd d (see right 
panel in Fig. 6), we see that we have h(i) = h{i + 1) for odd i, i.e. that the correlation function decays in steps of 



14 



two. This can be understood by the following argument. Consider a system of sequences s' = (s,s„+i) G H^^^, i.e. 
with the constraint that d'{s',0) is even, where d'{s\0) = d{s,0) + s„+i is the distance in the Hamming space Hl^'^^ ■ 
This constraint fixes the last bit Sn+i^ as s„+i = if d{s, 0) is even and Sn+i = 1 if d{s, 0) is odd. If we consider in 
the original space H2 a sphere at distance d{s, 0) = i with i odd we have d'[s' , 0) = i + 1, but also if d{s, 0) = i + 1 
we have d'(s', 0) = i + 1. Thus, spheres at distance i and i + 1 with i odd have the same distance i + 1 (even) from 
the origin in Hl^^^ . 

This means that the original problem on iJj with minimum distance d odd can be mapped into a problem on 
with d' = d+1 even. However, in the new problem, particles that were at a distance i and z + 1 from the (old) origin 
(for odd i) are at the same distance from, the (new) origin, thus the probability to have a particle at distance i (odd) 
and z + 1 in the original problem, given that there is a particle in the origin, should be the same. 

Note that the argument cannot be repeated for even d as in this case we should map the problem into a problem 
on H^^^ with distance d' odd, and this is inconsistent as discussed above. 

IV. CONCLUSIONS 

We discussed the high-density behavior of a system of hard spheres on the hypercube iJ^ from a physical point of 
view, trying to understand the mechanisms that determine the maximum density of the system. 

First we found a possible asymptotic solution for the liquid correlation in the limit n ^ 00, and we showed that 
this solution yields a negative entropy for ip <^ log(n)/n, i.e. for a number of particles N ~ |^n^/", very close to the 
Varshamov-Gilbert lower bound. 

On this ground we argue that a phase transition must exist towards a different phase. The nature of this transition 
is still unknown, but we presented some arguments in favor of a glass transition (basically, the analogy with the 
problem of hard spheres in M" for n ^ 00 [9]) and for a first order transition toward a phase in which the particles 
are constrained to a sublattice i?"g for even d (some insight from the solutions of the HNC equations and from 
preliminary Montecarlo simulations). 

Both these possibilities require further investigation. In particular, the study of the glass transition requires the 
computation of the replicated partition function of the system following [9, 15]. This is more difficult in this discrete 
problem as we cannot make a Gaussian ansatz for the single particle density. The study of the first order transition 
will require more extensive numerical simulations. 

Moreover, there is also the possibility that, at least for some particular values of n and d, some particular "crystal- 
like" configurations of high density exist. Unfortunately, our approach is based on a low density expansion (the virial 
series) so it is unable to capture the existence of such special configurations. 

The results presented here are not conclusive but in our opinion may lead to new ideas on how to improve the 
current bounds on R{6). We hope, in particular, that future work will clarify whether a glass transition exist or not 
in this system. 
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APPENDIX A: MORE ON SYMMETRIES 

We will discuss here the reason why we must restrict to even values of i in computing hi{x), see Eq. (55). As 
we discussed in section HIE for odd d the function h{i) has the property that h{i — 1) = h{i) for even i. This is a 
consequence of a hidden symmetry, namely the possibility to map the problem in iJ^ in a problem in H^~^^. 

It can be shown that if h{i) has such a symmetry, its Eourier transform h{a) has the symmetry h{a) — h{n~a + 1), 
a = 1, • • • ,n, due to the structure of the matrix J^„(a, i). This is consistent with the observation that in the limit 
n ^ 00 with a = nk, h{k) = h{l — k), as comes out from Eq. (58). However we are discarding a factor 1/n. 

When inverting the Fourier transform to recover h{i), we find that if i is even we get a meaningful result, while if i 
is odd, the symmetry h{k) = h{l — k), if interpreted as h{a) = h{n — a) leads simply to h{i) = due to the symmetry 
(18). The direct computation of h{i) for odd i would require a more refined calculation, however it is much simpler 
to compute h{i) for even i and use the identity h{i — 1) = h{i). 
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This procedure seems to produce meaningful results as evidenced by the positive agreement with numerical data, 
see Fig. 5. 
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